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We simulate self-reproducing micellar systems using a recently introduced lattice-gas 
automatonn. This dynamical model correctly describes the equilibrium and non-equilibrium 
properties of mixtures of oil, pwater and surfactants. The simulations reported here mimic 
the experiments of Luisi et ala in which caprylate micelles are formed by alkaline hydrolysis 
of immiscible ethyl caprylate ester. As in the laboratory experiments, we find an extended 
induction period during which the concentration of micelles remains small; thereafter the 
ester is consumed very rapidly with concomitant production of micelles. 

I. INTRODUCTION 

This article reports on simulations of autopoietic self-reproducing micellar systems that we have carried out 
using our recently introduced lattice-gas model for self-assembling amphiphilic systemsEJ. A self-reproducing 
system is one in which the process leading to population growth, here of amphiphilc and micelles, occurs as 
a result of the parent structures themselves. Additionally a structure which is self-bounded and is able to 
self-generate due to reactions which take place at or within that boundary, meets the criteria of autopoiesis, 
which certain authors have advocated as a minimal definition of lifeBu. 

Our model simulates the behavior of various experimental systems that have been investigated in recent 
years. Experiments on micelles that can catalyse their own reproduction were first described in and around 
1990-1991013. In these, micelles were present in the initial reaction mixture - that is, the system was 
presented with the bounded structures required for autocatalysis. Then in 1992 Bachmann, Luisi and LangEl 
reported on a two-phase system in which autocatalytic micelles were formed from amphiphiles that were 
themselves generated from a hydrolysis reaction within the aqueous solution. The system contained an 
aqueous sodium hydroxide solution and an ester, ethyl caprylate, which is itself immiscible with water. 
During the reaction, the system is well stirred in order to enhance the mixing of the two liquids. Hydrolysis 
of the ester initially produces the amphiphile monomer, sodium caprylate, which is known to form micelles 
in watem. This process occurs at a "slow" background reaction rate, largely confined to the ester-water 
interface. However, as soon as sufficient caprylate is present in the system for aggregation into micelles to 
take place (which occurs at and beyond the critical micelle concentration, in the terminology of equilibrium 
thermodynamics), a sudden increase in the reaction rate is observed. Bachmann et al. interpreted this as 
being due to a micellar catalytic process: When a micelle is formed, its hydrophobic interior causes it to 
rapidly bind oily ethyl caprylate molecules, dramatically enhancing the solubility of the previously immiscible 
liquid. The rate of hydrolysis is suddenly increased because the ester is now dispersed throughout the aqueous 
medium and its probability of encountering water-soluble hydroxide ions is vastly greater. Moreover, new 
amphiphile is produced by a chemical reaction that takes place at or within the domain boundary of the 
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parent micelles, thereby qualifying as an instance of autopoietic self-reproduction. This amphiphile is then 
able to form more micelles, and the hydrolysis reaction accelerates. The latter reaction process is the "fast" 
micellar catalysis step, the rate of micelle production increasing as the micelle concentration increases. The 
reaction is completed when there is no longer any ethyl caprylate left in the system. Upon lowering the pH 
of the system, Bachmann et al. found that the product micelles could be reversibly converted into vesicles 
formed from closed caprylate bilayer membranes. These authors argued that the properties of such systems 
may have relevance to chemical mechanisms associated with the origins of life on Eartho. 
In outline, the physicochemical steps occurring can be written down as rate processes: 

EC^C, 

nC^C n , (LI) 
C n + EC — ► C n + C. 

where EC denotes ethyl caprylate, C the caprylate anion and C n a caprylate micelle (n being the number 
of caprylate monomers per micelle, estimateda to be approximately 63). The third step here represents the 
micellar-catalysed hydrolysis of ester. 

Although the qualitative explanation of the Luisi experiments proffered here is simple enough to compre- 
hend, producing a quantitative model of the phenomenon is quite a different matter. The difficulties stem 
from the need to describe the dynamics of micelle formation correctly. Only recently has a detailed nonlin- 
ear dynamical model of-jpicelle formation kinetics been proposed and analysedEl, with particular reference to 
Luisi's system (see alscO). The approach is_hased on a generalisation of the Becker-Doring scheme describing 
the kinetics of first-order phase transitionsliil, in which micellar clusters grow (shrink) through the stepwise 
addition (removal) of individual surfactant monomers. The general kinetic process in the Becker-Doring 
theory can be written as 

a r ,b r+1 

C r + C C r +i, (1-2) 

where C r is a cluster comprising r monomers and r ranges from one to infinity. Here a r and b r +\ are the rate 
coefficients for the forward and reverse steps respectively. The reaction rates are given by the law of mass 
action, and lead to an infinite system of coupled nonlinear differential equations which may be combined 
with the processes shown in the scheme (O)B- 

At a more microscopic level, one would like to be able to describe micellar dynamics accurately. This is 
one aspect of the wider issue of the challenge of modelling self- assembling amphiphilic systems. Although 
in principle one might wish to use molecular dynamics to simulate the interactions in these systems in full 
molecular detail, in practice such an approach is computationally overbearing and will remain so for the 
foreseeable future given the state of today's technology. Moreover, some of the-jpost important situations 
demanding attention involve the coupling of self-assembly with hydrodynamicsE3. Their complexity is the 
result of the interplay between microscopic interactions and mesoscale hydrodynamic effects. Accordingly, 
we have developed a lattice-gas model of such systems which is capable of addressing these issues in a 
computationally tractable manner. 

The purpose of the present paper is to show how our model may be applied to simulate the dynamics 
of self-replicating micelles. Even at equilibrium, such structures are highly dynamic, rapidly exchanging 
monomers with the surrounding medium on timescales which are typically of the order of micro- to nano- 
seconds. This makes these processes very difficult to capture using molecular dynamics (MD), where the 
integration time-steps are on the order of femtoseconds. By contrast, lattice gas automata are based on a 
discrete time-stepping procedure performed on a spatially discrete lattice: The interactions are confined to 
collisions at the vertices, and the particles are then advected to neighboring sites in time steps that are on 
the order of the mean free time, several orders of magnitude larger than those in MD. 

The formation of micelles from surfactant monomers is an important part of the simulations described 
here. In addition, however, in the self-reproducing micelle experiments a chemical reaction is involved which 
converts the oily ester into surfactant molecules. Again, it is a strength of cellular automata that they 
enable chemical reactions to be included .easily by means of simple collisional transitioH-jrules; examples 
include the Belousov-Zhabotinski reactionllj, the Schlogl modeliil, and cement hydrationO-3. We shall find 
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that our suitably modified lattice-gas automaton model displays the same behaviour as that shown in the 
real system of self-reproducing micelles. It therefore bridges the divide between the two approaches to the 
investigation of autopoietic structures advocated by VarelaEil: virtual (computational) and physical (chemical 
synthesis) . Cellular automata were originally proposed as providing a virtual model of autopoiesis by Varela 
et aln; it is therefore fitting that we are able to connect the logical and the physical through a finite-state 
automaton simulation of the physical approach. Indeed, we hope that our results may be seen as being of 
wider significance than purely as simulations of phenomena in the natural world. For in recent years, there 
has been a renewed interest in the concept of "artificial life," the study of life as it could be (in whatever 
shape or form it may be found or made on Earth or elsewhere), rather than as it is known today on EarthEZI. 
According to this view, neither actual nor possible life is determined by the matter of which it is composed. 
For life is a process, and it is the form of this process, not the matter, that is the essence of life. As von 
Neumann originally sought-to demonstrate, one can ignore the physical medium and concentrate on the 
logic governing this processed. In principle, one can thus achieve the same logic in another clothing, totally 
distinct from the carbon-based form of life we know. 



II. DESCRIPTION OF THE MODEL 



Our lattice-gas model is a microscopic dynamical system which gives the correct mesoscopic and macro- 
scopic behaviour of mixtures of oiL water and surfactant. The model is based on the two-fluid immiscible 
lattice gas of Rothman and KellerEa, which we have reformulated using a microscopic particulate decription 
to permit the inclusion of amphiphile. The model exhibits the commonly formed equilibrium microemulsion 
phases, including droplets, bicontinua and lamellaad. Moreover, the model conserves momentum as well as 
the masses of the various species- and correctly simulates the fluid dynamical and scaling behaviour during 
self-assembly of these phasesElO. 

In order to incorporate the most general form of interaction energy (Hamiltonian) within our model 
system, we introduce a set of coupling constants a, [i, e, £, in terms of which the total interaction energy can 
be written aM 

AH mt = aAH cc + ^AH cd + eAH dc + (AH dd . {11.3) 

The four terms on the right hand side correspond, respectively, to the relative immiscibility of oil and water, 
the tendency of surrounding surfactant to bend around oil or water droplets, the propensity of surfactant 
dipoles to align across oil-water interfaces and the contribution from pairwise interactions between surfactant 
molecules. n 

In the experiments of Bachmann et a/.EI, the caprylate amphiphile produced has a tendency to form micelles 
in the bulk water phase rather .-than predominantly congregating in monolayers at the oil-water interface. 
Previous studies with our modelEJ enable us to simulate this property by choosing the following set of coupling 
constants: 

a = 1.0,/x=1.0,e = 3.0,C = 0.5, (II.4) 

for our simulations. There is another parameter present in our model that must be specified, called (3. This 
coefficient arises from the stochastic, Monte-Carlo-like, nature of the selection of outgoing collisional states 
at verticegj-within our lattice-gas automaton, and can be thought of as an effective inverse temperature-like 
parameters. 

From the discussion of Luisi's self-reproducing micelle experiments given in Sec. |, it is clear that the 
autocatalytic nature of the reaction is related to the critical micelle concentration (c.m.c.) of the caprylate 
anion. Experimentally, this is the point when the equilibrium concentration of surfactant monomers reaches 
a level at which they associate into micelles. However, it is important to recognise that the c.m.c. is not 
the location of a sharp phase transition, but rather the concentration at which the equilibrium fraction 
of monomers in micelles reaches some arbitrary value, usually taken to be 0.5. Indeed, the equilibrium 
fraction of monomers within micelles is itself somewhat arbitrary, since there is no clear-cut division between 
"micelles" and smaller (or, for that matter, larger) clusters that are not regarded as micelles. The equilibrium 
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fraction of monomers in micelles varies very rapidly with monomer concentration around the am. a, so 
that to an experimentalist it may look like a phase transition. As in our work based on the generalised 
Becker-Doring equationstJ, our lattice-gas model dispenses completely with the commonly made assumption 
of thermodynamic equilibrium between monomers and monodisperse micelles, being a highly dynamical, 
noncquilibrium model. Thus, we do not assume the existence of a critical micelle concentration: Rather, the 
phenomenon is itself an emergent property of our model. 

At the mesoscopic level, our model is very useful for analysing the kinetics of micellar growth. Indeed our 
model has already been shown to exhibit a critical micelle concentrat ion i n the binary - surfactant and water 
(or oil) - limitEJ. Using the set of coupling constants as defined in Eq. II. 4, we reproduce this analysis in order 



to obtain an estimate for the number of surfactant particles required to be in the binary water-surfactant 



system for the c.m.c. to be exceeded. The results can be found in Sec. |TJ. _ . 

Although an extension of our full microemulsion model to three dimensions is currently underwayEj, in 
this letter we shall be concerned only with the two-dimensional (2D) version. This is sufficient to capture 
the salient properties of the micellar self-production experiments; we expect the results in the 3D case to be 
qualitatively similar. 

In the basic simulation (as in the experiments) the system is initially comprised of two immiscible fluids: 
The majority phase we identify with the aqueous solution in the experimental set-up; the other, minority, 
phase we take to be representative of the ethyl caprylate ester. We then specify two chemical reaction mech- 
anisms designed to represent: (i) the "slow" pseudo-first-order alkaline hydrolysis of caprylate by hydroxide 
ions present within the aqueous medium producing surfactant, and (ii) the "fast" micellar-catalysed hydrol- 
ysis during which self-assembled micelles greatly enhance the solubility of ester within the water phase and 
thus lead to an increased alkaline hydrolysis. As in the actual experiments, the simulation reaches completion 
once all the ester in the system has been hydrolysed. 

Within our lattice-gas automaton, alkaline hydrolysis is described algorithmically by the conversion of 
an ester (oil) particle into an amphiphile particle if certain prescribed conditions are met; the system is 
otherwise allowed to evolve in the normal way (that is, with collisions conserving momenta and particle 
masses) and a record is kept of both the total number of surfactant particles in the system and the total 
number of micelles. (The precise manner in which "micelles" are defined in these simulations is specified 
below.) These data can then be used to analysG-,our simulations in a manner similar to that used for the 
experimental data generated by Bachmann et aln. 

In order to simulate the self-reproduction of micelles, we set up the model as follows: 

1. The 2D simulation domain with periodic boundary conditions in both dimensions is initialised with 
water and oil in a 6 : 4 ratio, the two bulk fluid regions being entirely phase-separated in accord with 
the experimental starting conditions. 

2. The inverse-temperature parameter (3 was set to the value 0.137. This value corresponds to the oil and 
water acting as immiscible fluids during the simulation (that is, below the so-called spinodal point at 
which the fluids become misciblc), but where the natural fluctuations within the system (which are 
larger near the spinodal point) encourage invasion of the individual bulk, phase-separated fluid regions 
by particles of the opposite kind. This has a similar effect to continuous stirring in the experimental 
system, since both enhance the ester hydrolysis rate, in essence by increasing the interfacial area. 

3. At each timestep and for every site x on the lattice which contains at least one oil particle, we perform 
the following analysis 

• "Slow" (pseudo-first-order) hydrolysis step: If five of the six nearest-neighbor sites are dominated 
by water particles, and less than two of these nearest-neighbor sites have amphiphile present (recall 
that there can be up to seven particles per site), then one of the oil particles on site x is converted 
into a surfactant particle with unit probability, prior to determining the outgoing state of the 
collision process from the look-up tables. This simulates ester hydrolysis in an overwhelmingly 
aqueous environment. 

• "Fast" micellar catalysis step: If five or more of the six nearest-neighbor sites contain one or more 
surfactant particles and five or more of the six next-nearest-neighbor sites have water particles 
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present, we define this structure as being equivalent to a micelle*. Within this bounded struc- 
ture a surfactant particle is created from an oil particle on site x with probability one, prior to 
determining the outgoing state of the collision process from the look-up tables. If there is more 
than one oil particle at site x then the oil particle actually converted is chosen at random; this 
prevents any bias from developing in the velocity given to the newly created surfactant particles. 
Chemically speaking, this step simulates the hydrolysis of an ester molecule solubilised within a 
micelle, by hydroxide ions in the surrounding aqueous environment. 

• Otherwise the automaton is updated according to the usual collision rulesQ. 

4. We keep track of the total number of surfactant particles present in the system as well as the number 
of micelles (as defined above) existing at each time step. Note that in addition to these "oil-in-water" 
type micelles, we also need to take into account similar nearest- and next-nearest-neighbor structures 
for which the central site is surfactant dominated. Otherwise the number of micelles present in the 
system would appear to drop to zero once all the oil in our system is converted to amphiphilc. 



III. SIMULATION RESULTS 



Initially we report on the analysis of the c.m.c. in a binary water-surfactant limit of our model. Simulations 
are on a 64 x 64 lattice with the initial condition being a random distribution of the water and surfactant 
particles; the density of surfactant in the system is altered in each case. When sufficient amphiphilc is 
present for the monomers to form micelles, these impart characteristic structure to the system which should 
be discernible in our simulations. To perform a quantitative analysis we keep track of the circularly-averaged 
structure functions S(k, t) of the surfactant densityo. The micelles themselves are dynamic structures and 
are not necessarily very long lived; in addition energetic considerations mean that individual micelles do 
not grow without limit, hence we do not expect to see evidence of micelles coalescing and growing in an 
unbounded manner. The results of our simulations are shown in Figs. |l| and |^. 

From Fig. [l] it is clear that for a reduced density of 0.01 surfactant (« 287 surfactant particles) there is 
no structure formation during the timescale of the run, so the amphiphile exists as monomers. In contrast, 
for 0.03 surfactant (« 860 surfactant particles) the figure shows structure formation early in the simulation 
which is maintained and does not appear to grow without limit, indicating that we have just reached the 
c.m.c. To further clarify the approximate position of the c.m.c, Fig. |^ shows the results for reduced densities 
of 0.02 (« 573 surfactant particles) and 0.04 (« 1147 surfactant particles) surfactant. The former of these is 
below and the latter above the c.m.c. for this system. These results indicate that the number of surfactant 
particles required to be in the system for the c.m.c. to be exceeded lies between 573 and 860. 



°*Note that in reality there are twelve next-nearest-neighbor sites on a hexagonal lattice. The six we specify here are 
those that lie along the directions defined by the six nearest-neighbor lattice directions from site x: This defines points 
around an essentially circular structure (which would be spherical in 3D), is computationally simple to implement 
and suffices for our present purposes. 
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s=0.01,time=0 
s=0.01,time=200 
s=0.01,time=1000 

s=0.03,time=0 
s=0.03,time=200 
s=0.03,time=1000 



FIG. 1. Temporal evolution of surfactant density structure function S(k, t) for binary water and surfactant mix- 
tures. The figure contains results for both 0.01 and 0.03 reduced density of surfactant. 



s=0.02,time=0 
s=0.02,time=200 
s=0.02,time=1000 

s=0.04,time=0 
s=0.04,time=200 
s=0.04,time=1000 




FIG. 2. Temporal evolution of surfactant density structure function S(k, t) for binary water and surfactant mix- 
tures. The figure contains results for both 0.02 and 0.04 reduced density of surfactant. 

The result of a typical simulation of self-reproducing micelles using the lattice-gas model described in Sec. |0] 
is shown in Fig. |^. The plot consists of the temporal evolution (in dimensionless lattice-gas automaton time 
steps) of the number of surfactant particles in the system as well as the number of micelles on a lattice 
of size 64 x 64. The comparison with the experimental results of Bachmann et al. is good: We see a 
relatively long initial induction period during which surfactant particles are produced by the "slow," pseudo- 
first-order hydrolysis step only; this is followed by an extremely rapid increase in the number of surfactant 
particles in the system which occurs at exactly the point at which micelles start to form. For, in accord with 
the estimate given above for the binary system, the amphiphile concentration at this point (« 600 — 700 
surfactant particles on the lattice) is at the critical micelle concentration. Figure |i] shows a sequence of 
snap-shots from the above described simulation at three different moments during the simulation: (a) shows 
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the situation shortly after the initial time, when fluctuations have caused particles from the immiscible fluids 
to interpenetrate; (b) depicts the mixture at the onset of the rapid phase in the reaction as amphiphile 
begins to be generated freely; and (c) displays the homogeneous end state of the reaction, where all the oil 
(ester) has been converted into amphiphile. 




FIG. 3. Amphiphile and micelle concentration as a function of time step in the lattice-gas automaton model. The 
triangles are for the micelles and correspond to the values indicated on the right hand vertical axis. 




FIG. 4. Three snapshots in time from the evolution of the lattice-gas automaton simulation. Ethyl caprylate 
molecules are black, water is red and anionic caprylate is green. 

Fig. |5| contains the result for an exactly analogous simulation to that described above but now the system 
is given a non-zero initial surfactant concentration. We use a reduced density of 0.005 which corresponds to 
about 150 surfactant particles randomly spread around the 64 x 64 lattice. The presence of these particles 
should induce an earlier onset of the rapid growth phase since the system should reach its critical micelle 
concentration more quickly. In accord with the experimental results of Bachmann et ala, this is exactly the 
behavior that we observe. 
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FIG. 5. Amphiphile and micelle concentration as a function of time step in the lattice-gas automaton model. The 
triangles are for the micelles and correspond to values on the right hand vertical axis. This case has surfactant in the 
system at time zero. 

Note that in Fig. ||, the calculation of the number of micelles present at each timestep is done with 
approximately 150 extra surfactant particles in the system than is the case for Fig. |[ This accounts for 
the perceived discrepancy in the number of micelles found in solution after the reaction has completed in 
each case. The system shown in Fig. |^ is found to have a slightly higher number of micelles because of the 
extra amphiphile monomers initially present in the solution. Again, by looking for the end of the "slow" 
hydrolysis step, we can estimate the critical micelle concentration for this system; as expected it is found 
to be approximately the same as in the previous simulation and so is also in agreement with the prediction 
from the binary water-surfactant mixture. 

It must be remembered that these lattice-gas simulations are of a microscopic, particulate nature, in con- 
tradistinction with our nonlinear dynamical modelsa. They are therefore subject to statistical fluctuations. 
Accordingly, we have also performed a series of simulations on progressively larger 2D lattices (64 x 64, 
128 x 128 and 256 x 256) in order to investigate the role of these fluctuations. At the same time, we have 
ensured that, in all cases, the simulations start from the same interfacial area to volume ratio for the oily 
ester (in 2D, this is more correctly stated as the interfacial length to area ratio), as this is recognised to 
be an important parameter in the experimental set-upo. (Experimentally, more rapid stirring leads to an 
enhanced ester/aqueous-phase interfacial area and hence a faster overall rate of reaction.) We carried out 
five such simulations for each size of lattice, differing only in the precise specification of the initial condition 
- the distribution of initial particle velocities within the bulk regions at timestep zero is randomised with a 
different seed for each run, leading to different fluctuation effects in each case. The results are summarised 
in Table I, which lists statistics for the time taken for the "slow" induction period to come to an end and so 
also marks the point at which the rapid reaction phase begins. 



TABLE I. Effect of lattice size on end of induction phase in lattice-gas automaton simulations of 
Luisi's experiments 



Lattice size 


Average timestep 
to of "fast" phase 


Standard 
deviation 


Furthest deviation 
from t 


64 x 64 


3827 


784.71 


1193 


128 x 128 


3821 


210.84 


302 


256 x 256 


3835 


144.52 


221 
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We find that the average timestep at which kick-off occurs is essentially constant across these simulations, 
while deviations from average behaviour are greatest for the smallest simulation box size and least for the 
largest. This is exactly what one would expect: As the number of particles in the system becomes larger, 
fluctuations are less significant and the average behavior (of an ensemble of similar systems) becomes ever 
closer to that observed experimentally. 



IV. CONCLUSIONS 

Using our recently introduced lattice-gas model which describes mixtures of surfactant, oil and waterEI, 
we have been able to simulate the experimentally observed kinetics of Luisi's self-reproducing micelleso. In 
the simulations, a characteristic induction period is followed, at the onset of the model's critical micelle 
concentration, by a very rapid production of surfactant and micelles. The induction period is shortened by 
the initial addition of surfactant, and by increasing temperature (which can also be taken to correspond 
to higher stirring rates of the fluids inside the reaction vessel) . The length of the induction period is also 
sensitive to fluctuations within the system; however, for ensembles of systems in which the interfacial area 
to volume ratio of the ester initially present is held constant, the ensemble average of the induction time is 
essentially independent of system size. 

Our simulations bring out very clearly one essential element in the self-reproduction of caprylate micelles 
during the aqueous alkaline hydrolysis of ethyl caprylate. This is the highly dynamic nature of the entire 
process, including the reversible self-assembly, and break-up of the micelles themselves, which are also central 
features of our previous Becker-Doring modeH. Moreover, "micelles" generally persist as clusters of a range-af 
sizes. The self-reproducing micelle experiment has been proposed as a paradigm of an autopoietic processEj, 
but it must be remembered that the micellar structures are themselves evanescent - they are continually 
breaking up and reforming. This is a property shared with living things, whose overall forms may remain 
largely invariant but whose detailed molecular constitutions are perpetually in a state of flux. 
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